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ABSTRACT 

We examine the outer Galactic HI disk for deviations from the 6 = 0° plane by constructing maps 
of disk surface density, mean height, and thickness. We find that the Galactic warp is well described 
by a vertical offset plus two Fourier modes of frequency 1 and 2, all of which grow with Galactocentric 
radius. Adding the m = 2 mode accounts for the large asymmetry between the northern and southern 
warps. We use a Morlet wavelet transform to investigate the spatial and frequency localization of 
higher frequency modes; these modes are often referred to as "scalloping." We find that the m — 10 
and 15 scalloping modes are well above the noise, but localized; this suggests that the scalloping does 
not pervade the whole disk, but only local regions. 

Subject headings: Galaxy: disk — Galaxy: structure — Galaxy: kinematics and dynamics — ISM: 
structure — radio lines: general 



1. INTRODUCTION 

Although the topography of the ga s disk of the Milk; 
Way has previo us ly bee n mapped tWgjjterjiautJ [L9JL, 
[ Henderson et all 1 19821: iBurton fc Te Lintel Hekker 
1986, and many others), its complete Fourier struc- 
ture has never been quant itatively described, though 
iBinnev fc Merrifieldl |l998) have approximated its 
three lowest frequency terms. Which spatial oscillation 
frequencies are most important? Is the scalloping a local 
or global effect? This paper constitutes an in depth 
analysis of the shape of the outer HI disk, as well as the 
first quantitative analysis of the scalloping. 

The large-scale warp in the gas disk of th e Milky Way 
has been known since 1957 llBurkel Il957t iKerrl Il957t 
lEHl |Kerr_ 



lWesterhoutHl957t iKerr et al.l ll957L The warp has a 
large amplitude, rising to a height greater than 4 kpc 
at a Galactocentric radius of 25 kpc in the northern 
data. It is also asymmetric; in the south the gas falls 
about 1 kpc below the plane before rising back to it. 
By observing other galaxies, it is possible to develop 
a gene ral understanding of how warps behave. iBosmal 
(1991) found that at minimum half of all galaxies are 
warped, and that galaxies with smaller da rk matt e r halo 
core radii are less likely to be warped. iBriggsl (1990) 
claimed that a warp's line of nodes starts out straight, 
and at a transition begins to advance in the direction 
of rotation, with some exceptions. Another survey has 
shown that warps are common in galaxies with HI disks 
that are extended compare d to their optical compone nts, 
and are often asymmetric l)Garcia-Ruiz et aT1l2002bj) . In 
the Milky Way, many components other than HI also 
participate in the warp A partial list includes dust 
(jFreudenreich et al.111991 CO dWouterloot et alJ fl990L 
solar nei ghborhood stars llDehnenll998D . and IRAS point 
sources IDiorgovski fc Sosinll1989lL ~ 

Many efforts have been directed toward understand- 
ing the warp on a theoretical basis. Bending modes 
have long been suspected as the mechanism creating and 
maintaining the warp. Early work studying the evolu- 
tion of bending mode oscillations was hampered by a lack 
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of knowledge regarding galactic halos, but showed that 
the shape of the density falloff near the edge of the disk 
plays an important role in the sta bility of bending modes 
i|Hunter fc Toomrd IT969t iToomrd fl983^ . The distribu- 
tion of matter in the halo controls the ability of the disk 
to sustain long-lived bending wave warps, so studying the 
prope rties of the warp will tell us a bout the shape of the 
halo iBinnevI Il978t iSnarkd 119841) . iSnarke fc Casertanol 
(1988) argued that bending mode oscillations are plau- 
sible when the halo an d the self-gravity of the disk are 
taken into account, but IBinnev et al.l l)1998j) showed that 
in this situation the warp will wind up within a few 
dynamical times. Previous work has largely been con- 
cerned with the existence and behavior of m — 1 warps 
only, though Sparke] l) 19951) demonstrated the stability of 
771 = modes in an axisymmetric halo. 

Several other mechanisms have been suggested as 
possibilities for creating an d maintaining a warp 
l)Kuiiken fc Garcia-Ruizl I2001D . Gravitational interac- 
tion with satellites such as the Magellanic Clouds is 
a promising candidate ijBurkd Il957t lWeinbergl Il998t 
Wein berg fc Blitzll2005t) . but there is a longstanding de- 
bate of whether tida l effects are strong enough to produce 
the observed effect ()Kerrlll95^ . Indeed, even in galaxies 
with companions similar to the LMC, tidal amplifica- 
tion may not be strong enough to a ccount for the size of 
the warp ijGarc ia-Ruiz ct al. 20023). Ac cretion of matter 
onto the halo is another plausible cause l)Jiang fc Binnevl 
1999), or matter can accr ete directly onto onto the disk 
and torque the gas orb its ijLonez-Corredoira et aUfeOOfl 
Sanchez-Salcedo 2004). Intergalactic magnetic fields can 
act on a slightly ion ized gas disk to produce a warp 
l)Battaner et al.lll99flL or the intergalactic medium c an 
excite a warp through a wind l)Kahn fc Woitierlll959ft . 

Throughout this paper, we will refer to shorter wave- 
leng th (5-25 kpc) os cillations in azimuth as scallop- 
ing. iGum et alJ l|196(1h first noticed a "waviness" in the 
gas layer of the inner Galaxy; observational evidence 
for scallo ping turns out to be p rominent in the oute r 
Galaxy llHenderson et all Il982t iKulkarni et"all 11982ft . 
iFleckl ()1983ft investigated the possibility that the Milky 
Way scalloping results from the Kelvin-Helmholz insta- 
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bility, and iSparkel l|1995j) suggested that the scalloping 
should only be present in the ou ter parts of Galactic 
disks. iSpicker fe Feitzingerl l)1986|) find azimuthal cor- 
rugations in HI and other components along the spiral 
arms after removing the wa r p. In an N-Body simula- 
tion, lEdelsohn fc Elmeereenl l(1997t) found evidence for 
spiral corrugations due to gravitational interaction with 
a satellite galaxy. In this paper we will only investigate 
oscillation in Galactocentric azimuth, and not in radius. 

In §2 we describe the method of transforming the He- 
liocentric data cube into Galactocentric coordinates, and 
present maps of the surface density, average height off the 
6 = 0° plane, and vertical thickness of the outer Galaxy. 
While these maps are constructed using new data, they 
are not significantly different from previous work. In §3 
we perform global and local analyses on these maps to 
better understand the warp and the scalloping. 



2. METHOD 



2.1 



Data Processing 

We use the 21 cm Lei d en/Argentine/Bonn ( LAB) 
data llKalberla et all 120051: iHartmann fc Burtonl E" 



iBaiaia et alTl2005t lArnal et all 1200(1) to conduct a 
quantitative study of the warp and scalloping. The 
LAB survey is a combina tion of the LDS data set 
l)Hartmann fc Burtonlll997|) with South ern sky observa- 
tions from the IAR l(Arnal et alJ 12000(1 : however, much 
attention has been paid toward ensuring a uniform data 
set. The data are corrected for stray radiation. The com- 
bined survey maps the entire sky within —450 < v r < 400 
km s _1 with a resolution of 1.3 km s , this velocity 
range includes all of the gas in the Galaxy in circular 
rotation. We used the Hanning smoothed data, which 
have a velocity resolution of 1.9 km s _1 , and used only 
data with |6| < 30°. The signal of the warp in t he north 
can b e weakly traced beyond this elevation limit l)Burtonl 
1985), but the vast majority of the warp signal is included 
within our b range. 

The LAB survey contains a large number of angu- 
larly small emission features, especially at high latitudes, 
which are not contiguous with the gas in the disk, and 
may not even be in circular rotation. Features like these 
are particularly troublesome at larger Galactic latitudes 
since they can contain enough gas to contaminate our cal- 
culations, especially at large Galactocentric radii. These 
objects are removed with a median filter so we can focus 
on the gas in the disk. Points with brightness tempera- 
ture Tb < are temporarily filled with the value 0.01 K 
for this filtering only. We then cycle through the LAB 
data cube and calculate the median of each point and its 
12 nearest neighbors in £ and b at the same line-of-sight 
velocity v r (a two dimensional diamond shaped filter); 
call this median T m (£, b, v r ). Any point with Tt > 10 T m 
is replaced with 10 T m . Two examples of objects caught 
by this filter are M31 and NGC 6822. 

The survey grid in £, b, and v r is not convenient for an- 
alyzing Galactic properties. Ideally, our measurements 
would be equally spaced in the Galactocentric cylindrical 
coordinates R, <f>, and z; we must interpolate a Galacto- 
centric grid from the LSR-centered data. The Galactic 
azimuth (f> is defined such that it converges with galac- 
tic longitude £ at large R. To convert from observed 
£, 6, and v r to i?, 4>, and z we use the following mapping 



functions: 

£ = sin" 1 

& = tan _1 



R 

z 



= sin £ cos b 



§eOR)-e 



Rn 



1/2 



+v n (R) cos (f> 1 1 - sin 2 £j cos b. ( 1 ) 

Here, 0(i?) is the Galactic rotatio n curve, which we as - 
sume to be 220 km s _1 everywhere ((Brand fc Blitzl l993). 
Oo and Rq are 220 km s _1 and 8.5 kpc, respectively, r is 
the vector connecting the Sun's location to the point un- 
der investigation; r' is the projection of this vector onto 
the plane of the disk (with magnitude r'). 

Equation ^ is the transformation that results from an 
assumption of elliptical gas orbits with major axis along 
the 4> = 90°, 270° line. To derive this, we assume the gas 
moves on the orbit: 



u = 22O km s _1 
v R = v n (R) cos0. 



(2) 



vr is the magnitude of the velocity in the Galactocentric 
radial direction, vjj is a parameterization of the clliptic- 
ity of the orbit, which is free to vary with Galactocentric 
radius; we discuss our method of calculating vu{R) in 
the Appendix. At all points, vr/v^ < 0.1. These equa- 
tions are simply the epicyclic approximation for an orbit 
with epicyclic frequency 1 and the angle of ellipse ori- 
entation fixed. Although the gas orbits in the Galaxy 
are not likely to correspond to the fixed ellipse orienta- 
tion we describe, this configuration minimizes t he cor- 
rection to v r for gas far from £ = 0° or 180° (see iVoskesI 
lfl99l for more detail on elliptical gas orbits). With- 
out a correction for radial motion of the gas, there is a 
large asymmetry between the surface densities at Galac- 
tic longitudes on eithe r side of £ = 0° and £ = 180° 
l(Henderson et alJ ll982). This must be taken into ac- 
count, or features in these two regions such as the surface 
density will appear discontinuous and distorted. Assum- 
ing an outward velocity for the Local Standard of Rest 
(LSR ) will correct this discontin uity to some degree 
1962; Kerr & Wcst erhoutl [l965'l . However, it seems that 
the best fit for the motion relative to the LSR changes 
with radiu s, implying that the effect is global, rather 
than local ((Blitz fc Spergelll991(l . To reduce the magni- 
tude of the discontinuity using gas orbits, one can use a 
Galactoce ntric radial velocity roughl y of the form cos <j> 
(or cost) l(Kuiiken fc Tremainell 19941) . 

We exclude all points that lie within 345° < £ < 15° 
or 165° < £ < 195°. Points in these two wedges have 
velocities along the line of sight that are too small with 
respect to their random velocities to establish reliable 
distances. All points in this region are set to Tt = 0. 

Using JU we construct a Galactocentric grid Tb(R, <p, z) 
by trilinear-interpolating from the grid Tf,(£,b,v r ). We 
do this by calculating the coordinates of an (R, 4>, z) point 
in {I, b, v r ) space, and interpolating from Tf>(£, 6, v r ). The 
resolution of the Galactocentric grid is set by the spacing 
of the LSR centered grid, but in this paper we are not 
interested in small-scale disk structure. A grid of 100 
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points in 10 kpc < R < 30 kpc, 350 points in — n < 
4> < it, and 141 points in -20 kpc< z < 20 kpc gives 
us sufficient resolution to answer the questions we are 
interested in. 

Undcrsamplcd grids fail to utilize all of the information 
in the data; our choice of grid spacing is both is an un- 
dersampling and an oversampling of the information in 
Tf,(£,b,v r ) depending on the position in the disk. Con- 
sider two points near R = 10 kpc and I — 15°, where 
we have poor Galactocentric resolution in (f>. The LAB 
survey has Ai = 0.5°; at this location this corresponds 
to A(/> w 0.9°. The spacing in our Galactocentric grid is 
larger than 1°, thus we are undersampled everywhere in 
the 4> dimension. In i)3.1.2l we will perform an azimuthal 
frequency analysis of each ring; none of the frequencies 
we examine approach the Nyquist frequency of the data. 
Near the midplane at R = 30 kpc and I = 15°, lines 
of constant b are separated by around 300 pc, and near 
the top of the grid they are separated by around 400 
pc; there is some oversampling in the z dimension by no 
more than a factor of 2. The most severe case of over- 
sampling is in the R dimension, where along £ = 15° we 
have only 14 velocity resolution elements in our R range 
for the 100 grid points. On the other hand, at I = 90°, 
there are more than 60 resolution elements in our range. 
We have chosen a grid spacing in R that oversamples to 
a varying degree depending on location in the disk. 

From Tb(R, (f>, z), we can recov er p(R,(f>,z) using the 
method outlined in iKerrl l)1968|) . We assume T s , the 
spin temperature, is 155 K everywhere. This is slightly 
higher than the maximum brightness temperature found 
in the LAB survey in the region we are concerned with. 
Choosing a larger number to force the optically thin limit 
(Tb <C T s ) makes a difference only in a small number of 
areas in the inner radii of our grid. The vast majority 
of the points are optically thin with any reasonable T s , 
and are not affected by this choice. The transformation 
to a density grid depends on \dv r /dr\. Using elliptical 
orbits makes calculating this quantity slightly more dif- 
ficult than with a flat rotation curve. Although the full 
derivative can be written analytically, we just calculate 
it numerically. Points with T\> < are set to p = 0. 



2.2. Surface Density and Mean Height maps 

The grid p{R, 4>, z) contains information about the den- 
sity of HI in the Galaxy, minus whatever has been re- 
moved by the median filter and the excluded regions. 
Previous studies have proceeded by calculating a mean 
height z(R, <f) for the gas. However, the Galaxy is a com- 
plicated place that contains a variety of HI structures in 
addition to the disk. In particular, there are many ex- 
tended clouds located near the disk as well as spurs that 
split off from the disk. None of these will have been re- 
moved by the median filter, whi ch acts only on compar- 
atively small areas of the sky. iVoskesI l|1999ft developed 
an alternative method to remove some of these features; 
he masked out a map of high velocity cloud complexes. 
Since we are only interested in the shape of the disk it- 
self, these additional components must be filtered out 
before our calculation of the mean height. We perform a 
dispersion filter that operates as follows. 



1. Calculate the total mass surface density 

JV=141 

M(R,<f>) = /W,*i)Az (3) 
where Az is the z bin size. 



2. Calculate 
z(R,c 



Sili 141 Zjp(R : 0, Zj)Az 
M(R,4>) 



3. Calculate the second moment 



M(R,cp) 



d 2 (RA) = 



(4) 



(5) 



4. Run through each point in (R, <p, z) space. For any 
point that does not lie within 2d of z, set p = 0. 
Call the grid that results from this dispersion filter 
p d (R,<p,z). 

For step 4, we experimented with several different cutoffs 
(in the range 1 — 3d) ; our results do not depend strongly 
on which cutoff we choose. 

We construct the Galactic disk surface density and 
mean height maps from pd{R,4>, z). For example, we 
can sum the density over the z dimension to construct 
the surface density associated with the disk (Figure QJ: 



E(R,<f>) = ^p d (R,<P,Zi)Az 



(G) 



The indices t and b represent the top and bottom Zi that 
have not been zeroed out by the dispersion filter. The 
resulting figure clearly demonstrates the falloff of the disk 
surface density with radius. The contour lines are nicely 
continuous across the I = 0° and t — 180° lines because 
of our use of elliptical orbits; see the Appendix for a 
version of this figure without correcting for these orbits. 
The jagged nature of the contours with R < 18 kpc is 
due in part to spiral arms l)Levine et al.ll2006|K There is 
a region near R « 27 kpc and a 90° with a smaller 
surface density than other regions at the same radius; 
this region has somewhat unusual features in all of our 
maps. Also notice the excluded regions near the Sun- 
Galactic center line; these gaps will appear in all of our 
plots. 

We will be looking closely at the mass weighted mean 
height of the gas disk (Figure [2J) , 



h(R,, 



1 



S(i?,. 



y^Zjp d (R, <j>,Zj)Az. 



(7) 



This height is calculated with respect to the Galactic 
midplane defined by b — 0°. The Galactic warp is the 
most immediately evident feature in this map; the gas 
in the northern hemisphere peaks at h i=a 5 kpc, while 
the southern gas descends only to h w —1.5 kpc, consis- 
tent with previous maps of the Galaxy. At least three 
vertical oscillations of magnitude w 1 kpc can be seen in 
the south from <j> w —120° to </> as —20°; these have pre- 
viously been called the "scalloping" . The approximate 
extent of the scalloping is marked with an arc connecting 
two "S" labels. The region of low surface density noted 
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in the discussion of Figure^(^ = 27 kpc, = 90°) nas 
a height that seems anomalous when compared to sur- 
rounding gas (it is marked with an "X"). Several features 
are elongated along lines of constant I indicating some 
level of contamination by turbulent velocities and/or lo- 
cal gas. In the past these have been dubbed "fingers of 
God" becau se they all point back to the Sun. 

in fro we will need a measure of the uncertainty in 
the mean height of each point in the map. We define this 
error using the sum of squared residuals, as is usual for 
least squares fits: 



E^fc'C* ~ h(R,<p)) 2 Pd (R, 0, z t )Az 



(8) 



(i-6 + l)S(i?,« 
where t — b + 1 is the number of points in the calculation 
of the mean and surface density. This is only an ap- 
proximation because it does not account for any of the 
uncertainty introduced in Tf,{£, b, v r ), vu(R), the inter- 
polation to p{R,<f),z), or the dispersion filter. Typically, 
e{R)/R^0M. 

Because of effects like turbulence, anomalous veloci- 
ties, and spiral arms, the features in the mean height map 
may not correspond to the actual shape of the Galaxy. 
The severity of this effect can be seen by looking at a 
contour map of dr/dv r (se e Figure |3"|) . While similar to 
previously published plots (Burton & Tc Lintel Hckkcrt 
1986), this figure also includes the correction for ellipti- 
cal gas orbits. Given a turbulent velocity of magnitude 
vt, features with coherent scales less than vtdr/dv r could 
potentially be false signals. Furthermore, even real fea- 
tures will be blurred out or even incorrectly positioned 
over the same length scale. A typical turbulent veloc- 
ity is around 8 km s . Small regions can differ from 
the fl at rotation curve by 20-30 km s _1 (Br and fc Blitzl 
1993). We do not expect that these distortions will sig- 
nificantly effect our Fourier analysis. The steep contours 
near £ — 0° and I = 180° demonstrate the need for the 
excluded regions, because small irregularities in velocity 
there result in large changes in distance. 

2.3. Thickness map 

We also construct a measure of the thickness of the 
disk. In contrast to the previous section, we keep all 
of the complicated features of the HI in our calculation 
of the thickness. We do this because these features do 
contain information regarding pressure and gravitational 
force, although they are a nuisance when calculating the 
mean height. In particular, our method for finding the 
thickness of the disk relies heavily on the tails of the 
vertical density distribution. Thus, in this section, we 
will work w ith p instead of pa . 

Following iHenderson et all l)1982j) , we define the first 
and third quartile points as Zji and Zjs as the smallest 
and largest indices, respectively, that satisfy: 
ji 

1 p{R,<j>,Zi)Az>M{R, 0)/4 



E< 

i=0 
141 

E 

8=j 3 



p(R,<p,z t )Az>M(R, (/>)/4. 



(9) 



We then define the half thic kness (a factor of 2 smaller 
than lHenderson et al.l l|1982f l to ease comparison with re- 
cent work): Th(R, 0) — (Zj3 — Zj\)/2. Note that this im- 
plies that Th is quantized by Az/2 (about 140 pc). This 



can lead to inaccuracies in Th in places where the thick- 
ness is small, i.e. R w 10 kpc. Th(R, 0) is shown in Fig- 
ure^ Some points near the sun actually have calculated 
half-thicknesses of zero; this is due to the comparatively 
small thickness of the disk in that region combined with 
our limited LAB survey range |6| < 30° and poor grid 
resolution in z. 

The flaring of the disk with radius is immediately 
evident, as Th increases in magnitude by a factor of 
about 8 bet ween R 10 and 30 kpc- this flaring was 
first seen in lLozinskava fe Kardashevl l)1963f) . Asymme- 
try between the northern and southern halves of the 
disk is also prominent; the southern half of the Galaxy 
has a lower average thickness t han the northern . This 
asy mmetry was also evident in [ Henderson et alJ l|1982f) 
and lBurton fc T c Lintel Hckkcrt (1986). The low surface 
density region around R k 27 kpc has a very large thick- 
ness; the gas in the region has clearly been disturbed and 
dispersed by some mechanism. There are regions of in- 
creased thickness close to £ = 15° and 345° near R w 10 
kpc probably caused by local gas. 

Figure [5] plots the surface density and half-thickness of 
the HI layer averaged over as a function of R. This 
figure confirms our visual impressions from the surface 
density and dispersion maps by showing the falloff in 
the surface density and the flaring of the thickness with 
radius. The contamination by local gas described in the 
previous paragraph as well as the poor z resolution prob- 
lem mentioned at the beginning of the section account 
for the rise in the average of Th with decreasing R near 
R « 10 kpc; there is no evidence the thickness of the disk 
actually behaves this way. A gra dient-expansion least- 
squares routine l)Marauardti 1 1 9631) gives the best fit for 
the surface density beyond 14 kpc (where the exponential 
falloff begins) as: 

S(i?) = 4.5xcxp[-(i?-14kpc)/4.3kpc] M© pc~ 2 . (10) 

3. ANALYSIS 

3.1. Global Structure 

Using the maps we have constructed, we conduct a 
quantitative investigation of the disk vertical structure. 
Since we are studying a disk, we use a method that is in- 
dependent of rotation in and treats = and <fi = 2n 
as the same point. Furthermore, the data are unevenly 
sampled due to the excluded regions within 15° of the 
Sun-Galactic center line. We complete the analysis with- 
out extrapolating h(R, 0) in these regions, to avoid in- 
troducing any artifacts into the signal. Throughout this 
paper, we refer to different frequency oscillations in the 
disk. These frequencies will always be labeled by the 
number of oscillations they will complete in a full 27r; 
thus the m — 1 mode has a one maximum and one min- 
imum in the disk. 

3.1.1. Galactic Warp 

The Galactic warp is the most prominent feature in 
Figure [21 A Lomb periodogram analysis of each radius 
ring (see q3.1.2f> reveals that the power in each of the 0, 1, 
and 2 modes is consistently larger than that in any other 
mode for R > 20 kpc. At some radii m = 3 is the next 
strongest mode, and at others it is m = 4. Accordingly, 
we characterize the warp by an offset in the z direction, 
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plus two Fourier modes with frequency 1 and 2. We fit 
each ring with the function: 

W(cf>) = W + Wi sin(0 - 0i ) + W 2 sin(20 - 2 ). (11) 

Each of the three amplitudes Wi and two phases <f>{ in 
this fit is a function of radius, because we fit each radius 
ring independently. We use the gradient-expansion fit- 
ting algorithm to perform this fit, weighting each point 
by the inverse of the squared estimate of the uncertainty 
in the mean height, e 2 (i?, </>). The results of this fit for 
the rings at R = 16, 22, and 28 kpc are shown in Figures 
16171 and|S| Error bars in these plots represent e(R,<f>). 



Each of these plots is a good fit; the offset and the two 
Fourier components are both necessary and sufficient to 
describe the large-scale structure of the disk. For illus- 
trative purposes, we will follow the ring at R — 28 kpc 
through each step of the analysis. 

Following the evolution of the fit parameters at differ- 
ent places in the disk will tell us how the warp changes 
with radius. The three amplitude parameters each in- 
crease monotonically, with the m = mode possibly 
reaching an asymptotic value near the far end of our 
radius range (Figure At R » 11 kpc, the m = 1 
mode dominates the shape of the warp; the other two 
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Fig. 2. — A contour map of h(R,<p). The shading changes every 500 pc in vertical displacement with respect to the 6 = 0° plane. A 
white contour is drawn at pc; black contours are drawn at 2000, 4000, and 6000 pc. The colorbar on the right hand side is marked in 
vertical kpc. The two "S" labels mark the approximate extent of the scalloping as seen by previous authors. The "X" marks the sharp dip 
discussed in A3. 21 and JJj 



modes do not become important until R « 15 kpc. This 
plot implies that the m = 1 mode has power even at the 
edge of our grid, thus we cannot establish the onset of 
the warp. We do a linear least-squares fit on the growth 
of each warp parameter using the function 

W n = k Q + ki (R - R k ) + k 2 (R- R k f (12) 

where only points at Rk and beyond are weighted in the 
fit. Rk is arbitrarily chosen to be near where each of the 
three modes starts growing. The value of ko for each fit 
is strongly correlated to the choice of Rk- 

iBinnev fc Merrifieldl l|1998t) (hereafter BM) discuss an 
approximation to the warp that is of similar form to our 
fit. They also fit the first three modes, but they fix the 
line of zeros to lie along the Sun-Galactic center line. In 
our fit, this would be equivalent to setting (f>i and <p2 to 
zero. Also, BM fix Wq and W2 to be the same. Note that 
our warp data are adjusted for elliptical gas orbits, and 



TABLE 1 

Parameters resulting from a linear least-squares fit 
to the warp 



m 


Rk (kpc) 


k (pc) 


fci (pc kpc 1 ) 


fe 2 (pc kpc 2 ) 





15 


-66 


150 


-0.47 


1 


10 


9 


197 


-3.1 


2 


15 


-70 


171 


-5.3 



are filtered according to the procedure described in 
whereas BM's data are not. Fig. ITTil compares the mode 
amplitudes calculated from the data, the fit from BM, 
and our fit. BM overestimate the strength of the warp 
starting at around R « 22 kpc for the m = 2 mode, 24 
kpc for the m = 0, and 27 kpc for the m = 1. However, 
the BM fit matches our data fairly well for the radii where 
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Fig. 3. — A contour map of dr/dv r , including the velocity 
contribution from the elliptical gas orbit correction. The axes are 
marked in kpc; the contours delineate regions of 250 pc/(km s _1 ). 
The contours rise close to the Galactic center and anticenter due 
to the small radial velocities of gas in those regions, demonstrating 
why the data in those directions are excluded from the rest of the 
study. 



the warp is growing most rapidly. 

The line of maxima for the m = 1 mode and one of 
the two lines of maxima for the m = 2 mode are marked 
on the contour map of the warp fit fFigure lTT|l . Since fa 
and fa are actually the line of zeros, the lines of maxima 
are shifted 90° and 45° from these values in our fit, re- 
spectively. The parameter fa is not well determined at 
small radii in our fit because the amplitude of the m = 2 
Fourier component is very small in that region. For this 
reason, we do not plot fa in the region where W2 is less 
than 150 pc. There is little evidence for precession in the 
lines of maxima for the two modes, and the line of max- 
ima for the m = 1 mode is roughly aligned with one of 
the lines of maxima of the m — 2 mo de; for each ra dius 
their difference in tf> is less than 12°. iBurtonl l)1988j) has 
examined how the two lines of zeros for the mean height 
change in <fr as a function of radius; they appear to stay 
roughly aligned with the Sun-Galactic center line as R 
changes from Rq to 26 kpc. Fig. II II shows that the line 
of zeros for the sum of our three warp modes falls within 
the excluded region, consistent with the earlier work. 

The three component fit shown in Figure 1111 does a 
good job of reproducing the large-scale features in the 
mean height map. This is not surprising, given that the 
rn = 0,1, and 2 modes are the strongest in the Lomb 
periodogram; differences between the two plots are due 
to power in high frequency modes. We will now examine 
the differences between the mean height map and the 
warp fit. 

3.1.2. Scalloping 

From the fit to the warp, we can determine a func- 
tion for the scalloping: s(R, fa — h(R, fa — W(R, fa. We 



continue to assume no information regarding the shape 
of the disk in the excluded regions. We can now look 
for higher m modes of oscillation that may be present 
in s(R,fa. Such a signal would be present if the gas 
were coherently moving up and down as a function of 
4>. The Lomb periodogram is a useful numerical method 
for detecting peri odic signals in unevenly spaced data 
ijPress et al.lll992T) . It provides the same amplitude we 
would determine by using a linear least squares fit, even if 
the signal is not sinusoidal, while allowing for a straight- 
forward error analysis. The data are weighted equally 
per point, which is necessary to deal with unevenly sam- 
pled data. To ease comparison with the warp component 
amplitudes in Fig. [51 we use an unnormalizcd Lomb pe- 
riodogram and take the square root of the power to get 
the amplitude. 

We use a Monte Carlo algorithm to determine the noise 
level for the Lomb periodogram. For each ring, the null 
hypothesis is Gaussian white noise with the same dis- 
persion as the data in s(R). We construct 10 3 sets of 
noise for each ring, run the Lomb periodogram on each 
set, and record the highest peak in the Lomb amplitude. 
We then determine the distribution of peaks in these am- 
plitudes, and define the 95% confidence interval as the 
amplitude just larger than 95% of the noise amplitude 
peaks. We use the same technique to determine the 99% 
confidence level. To conclude that a signal is real and 
not caused by noise, the signal strength must cross these 
confidence levels. We found nearly identical noise lev- 
els by scrambling the order of the data instead of using 
Gaussian white noise, indicating that these calculations 
are robust. 

The Lomb periodogram of s(R, fa for the R = 28 kpc 
ring is shown in Figure El 95% and 99% confidence 
levels are marked as thresholds in amplitude. At this 
radius, there is significant strength in modes 4-6, 10, and 
15. These modes are markedly weaker than the warp 
components at the same radii, which have amplitudes of 
1-2 kpc. 

In the case of unevenly spaced or missing data points, 
simple sine waves are not eigenmodes of the system. One 
way to see this is to take a Lomb periodogram of a 
pure sine wave on our <fi grid with missing data; a small 
amount of the power will leak into other frequencies. 
Thus subtracting out the warp from h(R, fa influences 
the amplitudes of the higher order modes from the Lomb 
periodogram because individual m modes are not inde- 
pendent. We argue that the subtraction is nonetheless 
acceptable because the warp and the scalloping appear to 
be due to physically distinct phenomena; studying them 
is much easier once they have been separated. Removing 
the strongest modes will also result in a more accurate 
power spectrum of the weaker frequencies, since we elim- 
inate the power leakage from the stronger to the weaker 
modes. 

Figure H31 shows how the amplitudes in several different 
m modes evolve with R. Modes 4-6, 10, and 15 appear 
to increase in strength at outer radii, while mode 3 is 
strongest in the intermediate radii in our map. Other 
modes have no detectable strength because they do not 
cross the 95% confidence level at any radii. Also notice 
how the confidence thresholds rise with radius because 
of the increase in the dispersion of s(R, fa. 

Another consequence of using the Lomb periodogram 
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Fig. 4. — A contour map of T^R, <f>); the thickness of the disk is approximately twice this quantity. Contours are drawn at 0, 500, 1000, 
1500, 2000, 3000, 4000, and 6000 pc. The colorbar and the figure labels are in pc. 



analysis on unevenly spaced data is that it is possible to 
be fooled by a false signal due to ringing from interfer- 
ence between different m modes. It is difficult to protect 
against this, but if two strong modes were interfering 
with each other to produce a third signal, we would ex- 
pect at least two of the signals to increase in strength 
at the same radius. The strongest modes in Figure 1131 
become significant at different radii and have visually 
different evolution with radius; we conclude that these 
modes are real. 

3.2. Local Structure 

The Lomb periodogram cannot be used to study the lo- 
cal structure of the disk since it cannot determine where 
in azimuth each mode is strong. Imagine a situation 
where, like a falling stone creating ripples on the surface 
of a pond, something passes through the HI disk and 
excites a local series of vertical oscillations. The oscil- 
lation will add power to some frequencies in the Lomb 



periodogram, but this effect may be dwarfed by oscil- 
lations elsewhere in the disk. We wish to detect these 
sorts of perturbations and learn where in the disk they 
are prominent. 

One way to draw out this type of structure is to use a 
wavelet analysis. Wavelets are ideal for our purposes be- 
cause they are localized in both position and frequency 
space; in other words, they will show where in the disk a 
given mode of oscillation is dominant. In wavelet theory, 
it is beneficial to choose a mother wavelet that is similar 
in shape to the signals you are searching for. Since we 
are looking for sinusoidal perturbations, we use the nor- 
malized Morlet wavelet with luq — 6; this is simply a sine 
wave times a Gaussian envelope. In the same way that 
a Fourier transform breaks down the frequency structure 
of a signal using sines and cosines, a wavelet transform 
breaks down a signal in terms of a set of these Morlet 
wavelets centered at different spatial positions and with 
a range of frequencies. However, in a continuous wavelet 
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Fig. 5. — The upper panel is the disk surface density, and the 
lower is the half thickness of the g function of radius. In the 

lower panel, the solid line is the average thickness over the whole 
disk, the dotted line is the average from 0° < <j> < 180° , and the 
dashed line is from 180° < </> < 360° . 
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Fig. 6.— The solid line is the warp fit W(<t>) to h at R = 16 
kpc. The approximate error in the mean height of the disk at each 
point is represented by the error bars. 



61 I I I 




-2 I I I I 

180 90 -90 -180 

Goloctocentric Azimuth (degrees) 

Fig. 7.— The solid line is the warp fit W(<f>) to h at R = 22 
kpc. The approximate error in the mean height of the disk at each 
point is represented by the error bars. 
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Fig. 8.— The solid line is the warp fit W{<j>) to h at R = 28 
kpc. The approximate error in the mean height of the disk at each 
point is represented by the error bars. 



transform like the one we perform, the different wavelet 
functions are not orthogonal. 

To avoid having to interpolate or zero out the ex- 
cluded regions, we will examine the northern and south- 
ern halves of the Galaxy separately. As in i Kt.1.21 we will 
work with the mean height function once the warp has 
been subtracted, s(R,<p). Points in frequency and posi- 
tion space that would involve the excluded regions are 
defined as being inside the "cone of influence"; points 
in this region are subject to edge effects, and are there- 
fore discarded. The points affected by this are not just 
those in the excluded regions but also those adjacent 
to these regions, within some range set by the wave- 
lcngt h (and thus the freq uency) of the particular mode 
l)Torrence fc Comr)dl998() . Thus, lower frequency modes 
will have a larger portion of the disk fall inside the cone 
of influence, and have to discard a larger range of points. 

It is also important to have some analytic measure 
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Fig. 9. — The evolution of the warp amplitude parameters as a 
function of radius. Wo is the solid line, W± is the dashed, and W2 
is the dotted. The dashed grey line is zero amplitude. 




Fig. 10. — In each pane l, th e solid line is the calculated ampli- 
tude of the the mode from 1111 . the dashed lin e is the fit from BM, 
and the dotted line is the fit described in 1121 . Top panel: m = 0. 
Middle panel: m = 1. Bottom panel: m = 2. 



of which peaks in the filtered power spectrum are sig- 
nificant. Significance lev els are discussed in detail in 
iTorrence fe; Compel l)1998|) . We construct a combination 
of parameters that have a x 2 distribution, and count 



as significant those that cross the 95% confidence level 
threshold. Again, we use Gaussian white noise to model 
the randomness in the height function. 

The wavelet power spectrum, W n , for the northern 
half of the R = 28 kpc ring is plotted in Figure the 
southern half is shown in Figure [TBI The mathematical 
details of the wavelet transform are summarized in the 
Appendix. This ring has several regions with significant 
power. These regions have a large width in both dimen- 
sions because wavelets do not have precise resolution in 
either position or frequency space. Much like the uncer- 
tainty principle, the cost of using a technique that gives 
both position and frequency information is mediocre res- 
olution in both. The power around <p ~ 90° comes from 
the sharp dip at that azimuth; as the figure shows, this 
causes ringing for a large range of frequencies. The power 
at -45° > > -90° is due to the "scalloping" in that 
region. 

The filtered wavelet power spectra for several different 
bands are plotted in Figure [TBI The band that is labeled 
m = 7 is actually the filtered sum of the power spectra 
that satisfy 6.5 <m< 7.5, and so on. Note that since we 
calculate the power spectra for the northern and south- 
ern halves independently, the significance contours are 
different for the two hemispheres, in addition to being a 
function of R. In practice, this occurs because the vari- 
ance of s(4>) for the northern half of the Galaxy is larger 
than that for the southern, and the strength required of 
the power spectrum to cross the 95% confidence level is 
directly proportional to the variance. 

4. DISCUSSION 

We grouped the m — 0,1 and 2 modes together in 
the warp because of their similar magnitudes in a Lomb 
periodogram. A close look at their dependence with R 
(Fig. [5J demonstrates that there are additional similar- 
ities. Both the and the 2 mode are near zero until 
about 15 kpc from the Galactic Center. The mode 
grows linearly from this point outwards while the 2 mode 
declines a small amount, and then grows approximately 
linearly by about 1.2 kpc. In comparison, the m = 1 
mode starts out fairly large (ss 300 pc), declines slightly 
and then grows by about an order of magnitude. Though 
the m = 2 mode does break off from this pattern at 
larger radii, all three components grow essentially mono- 
tonically, approximately linearly, and with similar slopes 
over the range 15 < R < 22 kpc. This may be a clue 
that their origins involve the same physics, and helps 
to justify our classification of higher frequency modes as 
scalloping. Studying the radial dependence of the warp 
inside the solar circle seems a fruitful way to learn more 
about these three components. 

The filtered wavelet power spectrum maps (Fig. Iltjfl 
are a representation of the scalloping in the outer Galaxy. 
Due to frequency-position uncertainty relations, it is im- 
possible to establish precisely the oscillation frequency of 
any local disturbance. This uncertainty manifests itself 
in the wavelet transform by perturbations that are only 
somewhat localized in frequency and position space, and 
therefore have some width in both. 

These maps demonstrate that the m « 10 scalloping 
found near £ ps 310° and R » 25 kpc mentioned in pre- 
vious work ijHenderson et alJll982t iKulkarni et"al1 ll982) 
is real. The strong power around m = 10 in the wavelet 
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Fig. 11. — Left panel: h(R, tj>). Right panel: The fit to the warp is plotted, along with the lines of maximum amplitude for the m = 1 
(dotted) and 2 (dashed) modes. These lines are marked W\ and W2 respectively. The white contour line denotes a height of pc; black 
lines mark the 2000, 4000, and 6000 pc elevations. 




Fig. 12. — The Lomb Periodogram for the ring at 28 kpc, once 
the warp has been removed. The 95% and 99% significance levels 
are drawn as dotted and dashed lines, respectively. 



transform is accompanied by a large amplitude in the 
Lomb periodogram for to = 10 at the corresponding 
radii. 

Wavelet transforms are also subject to the same 
problems that bedevil traditional Fourier transform ap- 
proaches. For example, there is a large amount of power 
around R w 30 kpc and w 90° in 8 < m < 15. This 
power is most likely not a result of scalloping in all of 
these different frequencies over some range of </>; instead, 
it is probably due to ringing. As in a Fourier decompo- 



sition, sharp changes in height will cause ringing in all 
frequencies of a wavelet transform. Indeed, the height 
map (Fig. EJ) does have an abrupt drop near cj> ~ 90° 
which could cause this ringing; we mark this region with 
an 'X'. The same feature can be seen after the warp has 
been subtracted in the ring at 28 kpc plotted in Figure 

m 

We use the local and global analyses in conjunction to 
determine whether there is any scalloping mode that is 
present over a full 2n ring in the outer Galaxy. This is 
important because it will determine whether the mecha- 
nism that causes scalloping operates on a global or local 
scale. For this task, both techniques are necessary be- 
cause even a pure to = 10 oscillation will have some 
width in frequency space when put through a wavelet 
transform. However, this same oscillation will have a 
sharp peak in the Lomb periodogram. Therefore, in order 
to state that some mode exists over a full 2-7T, we require 
both significant strength in the Lomb periodogram and 
significant power over a large range of <fi in the wavelet 
transform. This combination will conclusively determine 
whether the scalloping is a global or local phenomenon. 

For modes with to < 6 this technique is not useful. 
These modes are clearly important because their Lomb 
periodograms show significant power for 3 < to < 6. 
Unfortunately, these modes have large enough wave- 
lengths that the excluded regions always interfere with 
the wavelet transform regardless of where in the disk we 
look. A method of reducing the size of the excluded 
regions would alleviate this problem. However, this is 
likely to be difficult because eliminating the excluded re- 
gions would require a detailed knowledge of both vr and 
the distribution of gas velocities due to turbulence in the 



12 



Levine, Blitz, & Hcilcs 



m = 3 



m=4 



m = 5 




15 20 
R (kpc) 




15 20 25 30 
R (kpc) 



800 



— 600 



400 



200 




20 25 
R (kpc) 



m = 6 



m = 7 



m=8 




15 20 
R (kpc) 




15 20 25 30 
R (kpc) 




20 25 
R (kpc) 



m = 9 



m=10 



m=1 1 



800 



~ 600 



400 



E 200 




800 



~ 600 



400 



200 




15 20 25 30 
R (kpc) 



800 



~ 600 



400 



200 




20 25 
R (kpc) 



m=12 



m=13 



m=14 



800 



— 600 



400 



200 




800 



600 



400 



E 200 




15 20 25 30 
R (kpc) 




20 25 
R (kpc) 



m=15 



m=16 



m=17 




15 20 25 
R (kpc) 




15 20 25 30 
R (kpc) 




20 25 
R (kpc) 



Fig. 13. — The Lomb periodogram amplitudes as a function of radius, once the warp has been removed. The 95% and 99% significance 
levels are drawn as light grey and dark grey lines, respectively. 
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Fig. 14. — The filtered wavelet power spectrum for the northern 
half of the R = 28 kpc ring once the warp has been subtracted. The 
95% significance level is marked with a solid line, and the cone of 
influence is outlined with a dotted line. The function s(28 kpc, </>) 
is plotted in the lower panel. 



disk. Portions of these region are optically thick, which 
will make density determination impossible. 

With these constraints, the Lomb periodogram leads 
us to conclude that the modes to = 10 and 15 are the 
most fruitful places to look for scalloping that travels a 
full circle around the disk. Other modes do have inter- 
esting features, but their smaller Lomb amplitudes imply 
that we could be fooled by the imperfect frequency reso- 
lution of the wavelet power spectrum. For example, the 
modes 7 < to < 9 have significant wavelet power for large 
portions of the northern half of the galaxy, but none of 
these modes has significant signal for the corresponding 
radius range in the Lomb periodogram. 

For the to — 10 mode the only large region with sig- 
nificant power in the south is the one around I w 310° 
that we noted previously. The power in this perturbation 
has fallen below the 95% confidence level by I « 270°, 
implying that this scalloping is a local effect. 

The same appears to be true for the to = 15 oscil- 
lations, although this mode does have a region of high 
significance near R « 30 kpc over a large range of <j>. 
However, the northern part of this signal is the region 
where we believe ringing to play an important role. It is 
therefore possible that the to = 15 mode does have sig- 
nificant wavelet power over the full 2tt near R « 30 kpc, 
but the evidence is not conclusive because the northern 
part of the wavelet power is most likely not due to scal- 
loping. 




-90 -180 

Galactocentric Azimuth (degrees) 
Fig. 15. — The filtered wavelet power spectrum for the southern 
half of the R = 28 kpc ring once the warp has been subtracted. The 
95% significance level is marked with a solid line, and the cone of 
influence is outlined with a dotted line. The function s(28 kpc, </>) 
is plotted in the lower panel. 



With the exception of the to = 15 mode, no frequency 
we examine with the wavelet transform carries significant 
power around an entire ring along with a correspondingly 
significant Lomb periodogram strength. For this reason, 
we conclude that scalloping generally appears to be a 
local phenomenon. The modes 3 < m < 6 remain a 
possible exception, since we were not able to study them 
with a wavelet transform. 

It remains unclear what mechanism acts as the energy 
source for the scalloping. One possible cause is a mas- 
sive object passing through the disk that excites local 
vertical oscillations in the HI gas. Azimuthally traveling 
wavefronts can be created by the magnetic field threaded 
through the disk; if the field is primarily azimuthal in na- 
ture, vertical oscillations will have larger phase and group 
velocities in the azimuthal direction than in the radial di- 
rection. This could lead to scalloping such as that seen 
towards I w 310°. 

5. CONCLUSIONS 

We fit the global shape of the warp on our grid of 
concentric rings using a vertical offset and two sinusoidal 
modes. Outside of R « 20 kpc, each of these three modes 
has more power than any of the higher frequencies we 
look at. The amplitude increases with radius over our 
entire radius range for the 1 mode, and starting from 
around 15 kpc for the and 2 modes. The growth of the 
and 2 modes results in asymmetry in the warp; this 
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Fig. 16. — The filtered wavelet power spectrum for the outer Galaxy. Darker contours represent higher power, and the 95% significance 
level is marked with a solid line. The labels on the rings give the distance from the Galactic center in kpc. 



growth begins near where the stellar disk ends. The line 
of maxima of the m = 1 mode is essentially coincident 
with one of the lines of maxima of the m = 2 mode. 
There is little evidence for precession or winding of these 
two modes. 

A global analysis with the Lomb periodogram shows 
that each m mode evolves differently with radius. The 
most interesting include m = 3 — 6, 10, and 15, each of 
which start at small radii below the 95% significance level 
and then cross it further out in the disk. An analysis 
combining the global Lomb periodogram and the local 
wavelet transform shows that none of the modes with 
7 < m < 15 have strength over a full ring of the disk. 



Using a wavelet transform, we show that the scalloping 
observed by previous authors near I ~ 310° is real. We 
therefore conclude that the scalloping is a local effect. 
Lower frequency modes proved impossible to study with 
wavelets due to the presence of the excluded regions. 

Wavelet software was provided by C. Tor- 

rence and G. Compo and is available at URL: 
htt p : / / paos .Colorado . edu /r esearch / wavelets / 

(Torrence & Comoo 1998). Many thanks to Peter 

Kalberla for providing a copy of the LAB data set. 
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APPENDIX 

MEASURING RADIAL MOTION OF THE GAS 

This appendix describes our method for finding a functional form for vu (R) , the magnitude of the elliptical corrections 
to circular rotation. A naive surface density plot using only circular rotation has a large degree of asymmetry across the 
lines £ — 0° and 180° (Kulkarni et al.lll982(l . Figure IA17I shows the calculated surface density without any correction 
for elliptical gas orbits, but with the filters we describe in [21 Observational experience tells us that the surface density 
should not have discontinuities in cardinal directions; the elliptical orbits we describe in ^\ provide the strongest 
corrections in the directions where the cos </> term in is large, and smaller corrections elsewhere. The algorithm is 
based on matching the surface density of the HI disk on either side of the excluded region centered at £ — 180°. We 
draw confidence from the fact that this fit also does a good job matching the contours around I = 0° even though they 
are not included in the fit. 

In order to calculate a surface density, we must first know vu(R) because it enters into the derivative \dv r /dr\. Thus, 
determining vu(R) from anything connected to p is a circular problem. We must assume some form of vn(R) and 
check to see if it results in surface densities that are well matched. We choose the functional form: 

v n {R) = a +P -2 , (Al) 

and do not include a zeroth order term to ensure that vn(R) passes through zero for the solar circle. 
iRadhakrishnan &: Sarmal l)1980j) showed that the LSR does not have a radial velocity with respect to the Galactic 
center. 

Our algorithm follows these steps: 

1. Construct vn(R) for some combination of a and /3 

2. Interpolate from the LAB survey to find Tb(R,£, z) for a grid in Ri and Zj over the survey range [6| < 30° and 
155° <£< 165° or 195° < £ < 205° using v n (R) 

3. Sum over z to find the surface density S(i?, £) 

4. Average over £ for the two subsets of £ to find Ei6s(i?) and £195 (i?) 

5. Calculate a modified \ 2 statistic to determine how well the surface densities are matched 
The modified x 2 is defined as 



x 2 = E 



^165 (Ri) < (Ri) 



where the < operator returns the smaller of its two operands. 

Using this algorithm we search for the values of a and (3 that minimize the modified x 2 - We find good matches for 
a = 8.67 and (3 = -1.08. 

We also tried a weighting that only fit points beyond the perturbations of the spiral arms, i.e. R > 2Rq. Although 
this did change vn(R) by 25% or so, it had no qualitative effect on our other results. 

WAVELETS 

The continuous wavelet transform for a discrete series of points z n is given by 

A 1 r (n'-n)A(j)' 



W n (s) = Zn'lp 
ra'=0 



(Bl) 



where N is the number of points in the series, n is a position index, A4> is the spacing of the points in <f> space, and 
ip* is the complex conjugate of the normalized Morlet wavelet ijTorrence fc Compd rl998): 

$ (r,) = ^\ V2 tt- 1 ^ eW'e-" 2 / 2 . (B2) 

The scale of the transform is s; for the Morlet wavelet this is simply related to the wavelength by A = 1.03s. Evenly 
spaced points are necessary in order to use this transform, but since we are treating the two halves of the Galaxy 
separately, s(R, 4>) is split into two halves, each with equally spaced points. We perform this transform for the dense 
set of frequencies given by 

Sj =s y A i, j =0,1,..., J. (B3) 

Here, sq is the smallest scale that can be sampled, 2A</>. Aj is a measure of how densely we sample in scale space; 
because computation time is not large, we choose relatively dense sampling throughout: Aj = 0.0125. 
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Fig. A17. — A contour map of £(i?, </>) without any correction for elliptical orbits. Contours are evenly spaced in logio space, at 
(3, 1,0.3,0.1) M /pc 2 . 



From IjBll) . wc construct the wavelet power spectrum, |W„(s)| . We calculate the wavelet power spectrum for the 
dense set of frequencies, and then filter over a range of scales to find the power in a frequency band. The filtered power 
spectrum is given by 



_ 2 _ AjAt A \W n ( Sj )\ 2 
n " c s h 7* 



(B4) 



where Cs is a reconstruction factor that depends on the choice of wavelet. For the Morlet wavelet with uiq = 6. 
C s = 0.776. 
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